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SMOOTHED CORNERS AND SCATTERED WAVES 


CHARLES L. EPSTEIN* AND MICHAEL O’NEILt 


Abstract. We introduce an arbitrary order, computationally efficient method to smooth corners 
on curves in the plane, as well as edges and vertices on surfaces in The method is local, only 
modifying the original surface in a neighborhood of the geometric singularity, and preserves desirable 
features like convexity and symmetry. The smoothness of the final surface is an explicit parameter 
in the method, and the bandlimit of the smoothed surface is proportional to its smoothness. Several 
numerical examples are provided in the context of acoustic scattering. In particular, we compare 
scattered fields from smoothed geometries in two dimensions with those from polygonal domains. 
We observe that significant reductions in computational cost can be obtained if merely approximate 
solutions are desired in the near- or far-held. Provided that the smoothing is sub-wavelength, the 
error of the scattered held is proportional to the size of the geometry that is modihed. 
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1. Introduction. In the numerical solution of boundary value problems for par¬ 
tial differential equations an especially difficult case arises when the boundary of the 
domain has corners (in two dimensions) or edges and vertices (in three dimensions). 
Several groups have devoted resources to solving this problem and have made seri¬ 
ous inroads towards addressing these issues in the context of the classical integral 
equations of mathematical physics (acoustic and electromagnetic scattering, elastic¬ 
ity, etc.) [10, 7, 13, 37, 34, 33, 16, 11, 45, 51]. The resulting numerical schemes 
often involve the use of specially designed quadratures which handle not only sin¬ 
gular or weakly-singular integrals but also singular layer potential densities. These 
methods are based on several standard ideas in modern numerical analysis, namely 
low-rank approximations [19, 6], generalized Gaussian quadratures and adaptive re¬ 
finement [54, 11], and (semi-) analytic product integration formulae [36, 35, 32]. More 
recently, explicit exact forms of the solutions to actual layer potential densities were 
derived in [51]. All the numerical tools just mentioned are now well-developed and 
require minimal sophistication to use, but can still be time-consuming, or too special- 
purpose to implement. An approach that has not been investigated thoroughly at 
this time in the literature is that of solving an analogous scattering problem from a 
smoother geometry that is close to the original one. In what follows, by close we mean 
different only in a (small) neighborhood of the geometric singularities (e.g. corners in 
two dimensions). 
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Lacking, up to this time, is a reliable, systematic, computationally simple method for 
smoothing such irregularities that also retains desirable geometric features, such as 
convexity or local symmetries. In this note we discuss several methods for doing this, 
including a simple convolution method as well as introduce a new geometric method, 
particularly useful for regularizing surfaces in three dimensions. Our methods are 
tailored for use with polygons in two dimensions and polyhedra in three dimensions. 
However, because the modifications are done locally, this approach can be applied 
to more general shapes (namely curves which intersect at their endpoints) through 
composition with diffeomorphisms. Indeed, our method already employs such com¬ 
positions in the three dimensional case. 

Corner and edge rounding methods are useful for two reasons. First, in the context 
of the solution of scattering problems via integral equations, smoothing geometric 
singularities on a sub-wavelength scale provides a means by which to apply stan¬ 
dard numerical quadratures [1, 43, 34] for weakly-singular integrals along smooth 
boundaries, instead of the more complicated schemes required in the neighborhood of 
corners. In two dimensions, our numerical examples show that convergence is roughly 
first-order in the scattered field, both in the near- and far-fields. In three dimensions, 
reducing the number of discretization nodes is particularly useful because of the rel¬ 
ative cost of even the fastest solvers. State of the art, high-order accurate solvers in 
three dimensions include those by Bremer, Gillman, Gimbutas, and Martinsson [9, 8] 
and Bruno [14]. 

Second, since the schemes to be presented only change the geometry locally, they may 
lead to a new class of algorithms which can be incorporated into modern computer- 
aided design (GAD) and engineering (GAE) software packages. The regularity of the 
smoothed surface can be precisely controlled in the neighborhood of the singularity. 
Applications in fine-grained polishing of machined mechanical parts are straightfor¬ 
ward. This paper investigates the advantages and disadvantages of solving a scattering 
problem from a nearby smoothed geometry instead of the original non-smooth one. 

We organize the remainder of the paper as follows. Section 2 reviews standard in¬ 
tegral equation formulations of acoustic scattering phenomena, as well as both the 
analytic regularity results of scattering from geometries with corners and the numer¬ 
ical techniques that have been developed to compute them. Section 3 describes a 
straightforward and systematic way to smooth the corners of polygons in two dimen¬ 
sions. The method can be extended to regions with piecewise smooth boundaries 
via the application of a diffeomorphism. Several numerical experiments are presented 
to illustrate the heuristics of the approach. Section 4 details several methods for 
smoothing polyhedra in three dimensions. The methods of Section 3 are extended to 
three dimensions and a new geometric method is introduced which is applicable in 
most cases. Section 5 puts all the previous sections together and gives a recipe for 
smoothing a general polyhedron in three dimensions. Section 6 reviews some ana¬ 
lytical methods that can be used to construct the diffeomorphisms required by the 
three-dimensional methods of Section 4. Lastly, the conclusions in Section 7 discuss 
drawbacks and difficulties with our method, as well as points to future areas of re¬ 
search and applications. Numerical experiments are included throughout the paper 
to demonstrate the application of scattering from smoothed geometries, as well as to 
visually describe the results of the smoothing techniques. 
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2. Scattering in singular geometries. There are two questions that require 
answers when studying scattering (acoustic, electrostatic and electromagnetic, etc.) in 
singular geometries using integral equations. First, in the neighborhood of a corner or 
edge, what regularity can we expect in the solution for data with a given smoothness? 
And second, if a solution exists, which can he represented in terms of a layer-potential 
density, is the density continuous and how can it he numerically calculated? The first 
question has been studied in detail by Dauge, etc [23, 25]. The latter question is 
mainly an exercise in numerical integration, and has been thoroughly studied by Bre¬ 
mer, Bruno, Helsing, etc. See [10, 16, 37] for more details. Often, the numerical 
solution is a combination of sophisticated quadrature schemes coupled with an adap¬ 
tive discretization of the geometry (in order to correctly resolve complicated layer 
potential densities). We now give a very brief review of some results in both of these 
areas. 


2.1. An integral equation approach. Almost all of the classical partial differ¬ 
ential equations of mathematical physics can be reformulated in an equivalent integral 
equation form [31]. The integral equation form has many advantages, namely the di¬ 
rect handling of unbounded domains in the case where the solution of a PDF reduces 
to a boundary integral equation. Furthermore, when the integral equation is Fred¬ 
holm of the second kind, as is often the case, provable bounds exist on the accuracy of 
the solution which are directly related to the order of the quadrature rule used in the 
discretization [3, 2]. In this section, we summarize a basic Nystrom-type discretiza¬ 
tion of an integral equation for the Helmholtz equation that can be used to solve an 
exterior acoustic scattering problem. 

Time-harmonic acoustic wave propagation in homogeneous free-space (we address the 
two-dimensional version here) is governed by the Helmholtz equation, 

(2.1) (A + u{x) =0 in Mf, 


where u is related to the acoustic pressure and k is related to the wavenumber of the 
field, namely k = uj/c, where uj is the angular velocity and c is the speed of sound 
in the medium. In particular, often one is interested in the solution to a scattering 
problem in the presence of some inclusion 11, where the total pressure field is 
the sum of an incoming field and a scattered field u. If the boundary of the 
inclusion is given by F, then sound-hard scattering phenomena can be formulated as 
the following boundary value problem: 


( 2 . 2 ) 


{A ^ k^) u^^\x) = 0 


du^^^{x) 

dn 


= 0 


in \ H, 
on F, 


where d/dn represents the derivative with respect to the outward normal to F. This 
boundary value problem is also known as the Neumann scattering problem. Dirichlet 
boundary conditions = 0 along F correspond to sound-soft scattering problems. 
The solution to (2.2) is unique under a suitable decay condition, known as a Sommer- 
feld radiation condition, on the scattered field u. In particular, in two dimensions, 
the scattered field u must satisfy: 



u{x) — iku{x) 


= 0 , 


(2.3) 
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and djdr is understood to be differentiation in the radial direction. It is well-known 
that the Green’s function for (2.1) is given in terms of the zeroth order Hankel function 
of the first kind, and is normalized as: 

(2.4) g^{x) = \H'i\k\x\). 

Using this Green’s function, a solution to (2.2) can be expressed in terms of a single¬ 
layer potential 


(2.5) 


u{x) = Sk[cr]{x) 



y\)cr{y) ds{y), 


where s is arclength along F. After taking the proper limit as x ^ F from the exterior, 
this representation results in the second-kind integral equation for the density a\ 

(2.6) Icr + [cr] = - on T, 

or more explicitly. 


(2.7) 



" d 

drix 


9k{x,y) 


a{y) ds{y) = - 




for cc G F. 


The operator represents the normal derivative of a single-layer potential. If F 
is C^, then the integral in (2.7) is weakly-singular and can be evaluated using spe¬ 
cially designed quadrature rules [32]. There are several approaches to discretizing the 
continuous integral equation (2.7), namely Galerkin, collocation, qualocation, and 
Nystrom discretizations [24, 5]. The methods of this paper apply to all of these 
approaches (under suitable small changes); we briefly describe the Nystrom method 
for its simplicity. 

The Nystrom discretization of (2.7) replaces continuous functions and integrals by 
samples and sums of samples. Namely, for a given quadrature rule consisting of nodes 
and weights {xj^Wji} for the integral appearing in (2.7), we approximate the solution 
(T{xj) « CTj at each node Xj as the solution to the system of equations: 

1 g 

(2.8) —cTj + ^^Wj£ —— gk{xj^X£)(j£ = —u 

p 


for all j. Here we have explicitly shown that the quadrature weights can be a function 
of the outgoing node Xj . As the number of discretization points or order of quadrature 
increase, cij approaches the exact solution cF(xj). The previous linear system can be 
solved directly if the resulting linear system is small enough, or for larger systems 
using iterative (fast multipole methods and GMRES, etc.) [18, 50] or fast direct solvers 
[29, 38]. 

It should be noted that integral equation (2.7) fails to be uniquely solvable at a dis¬ 
crete set of /c’s, known as spurious resonances. This is not a failure of the uniqueness 
properties of the PDE, but rather a failure in the particular choice of integral repre¬ 
sentation. Ghoosing what is known as a combined-field representation can result in a 
uniquely solvable integral equation, albeit at the cost of a slightly more complicated 
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formulation [21, 14]. One possible combined-field (or regularized) representation of 
this type is of the form: 


(2.9) 


U = 5fe[cr] + ar>fe5o[cr], 


where a is a user-chosen complex-valued parameter, V is known as the double-layer 
potential, given by 


( 2 . 10 ) 





gk{x,y) 


a{y) ds{y), 


and iSo is a single-layer potential corresponding to the Green’s function for Laplace’s 
equation: 


(2.11) 5oM(a:) = ^ Tin^_-^cr(y)ds(y), 

There are many other regularizations that one may use, and this is the subject of 
ongoing research (especially in the large-/c regime). We make a point to explicitly 
state the form of the integral representation for numerical experiments appearing 
later in the paper. 


2.2. Analytic results in singular geometries. In the previous section we 
discussed the process by which the Helmholtz boundary value problem (2.2) for the 
field u is reformulated as a boundary integral equation for a separate unknown layer 
potential density a. We have not, however, discussed the effect that the geometry 
has on the solution a (assuming that the data is smooth). The regularity of the 
solution a to the integral equation is strongly affected by the presence of corners on 
the boundary T, the boundary data, and details of the local geometry, e.g. whether 
the corners are re-entrant, acute, obtuse, etc. 

On smooth domains, the layer potential operators 5/^, 5^, and are compact, 
classical pseudodifferential operators and therefore the invertibility of the associated 
second-kind integral equation follows from the Fredholm alternative [27, 28]. The 
mapping properties on Sobolev and Holder spaces are well-known and essentially 
optimal. However, when the domain is merely continuous and not everywhere dif¬ 
ferentiable, these operators cease to be compact. While canonical PDF results have 
existed for some time, it is a relatively recent result in functional analysis that the clas¬ 
sical integral equation corresponding to the interior Dirichlet problem for Laplace’s 
equation, namely 


( 2 . 12 ) 


:P{x) 


L 


d 1 


driy 271 \x — y\ 


In ■ 


p{y) ds{y) = f{x), 


for cc G r. 


where T bounds some Lipschitz domain D, is invertible on £2 [53]. Similar results 
exist for the Neumann problem as well, and, with some work, extend to the analogous 
integral equations in the Helmholtz case [22, 39]. 


Classically, representations for solutions to the Helmholtz equation can be obtained 
in the exterior of a wedge or corner by using fractional Bessel function expansions, 
as in [42]. An expansion of this type, however, does not immediately yield similar 
statements concerning the density, a. Very recently, however, expansions of the actual 
density (at least in the Laplace case) were derived that allow for the construction of 
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very efficient, most likely optimal, solvers [51]. The topic has been further studied 
by many in the finite element, asymptotics, and analysis communities including, but 
certainly not limited to, Buffa, Ciarlet, Costabel, Dauge, and others [25, 23, 17]. This 
classical work addresses solutions to the Helmholtz equation and Maxwell’s equations, 
as well as Stokes flow in fluid dynamics and elasticity. 


2.3. Numerical methods for Lipschitz domains. While the results of the 
previous section are interesting from a mathematical standpoint, and certainly offer 
insights on how to properly construct finite element methods that have desirable 
properties in singular geometries, they offer little help in the construction of numerical 
quadrature schemes that can be used efficiently in the Nystrom method solution for 
the associated boundary integral equation. Recently there have been several papers 
addressing the question of constructing (mostly brute force) discretization schemes 
for boundary integral equations on polyhedral domains or domains with corners. As 
mentioned before, these schemes are often a combination of adaptive refinement of 
the geometry near the singular set, the design of specialized quadratures, and proper 
re-weighting of the unknown density. 

Adaptive or dyadic refinement of the geometries and density near geometric singular¬ 
ities has been commonplace for some time, but it was only recently detailed how to 
embed the Nystrom discretization into the proper continuous function space in order 
for the spectrum of the finite-dimensional approximation to converge to the spectrum 
of the continuous integral equation [7]. We omit a discussion of the dyadic refinement 
methods since they are well-known and [37] offers a nice review. However, we briefly 
mention the £2 norm-preserving scheme discussed by Bremer. 

First, it should be pointed out that the unknowns in the discrete system (2.8) are 
point values of the continuous density a. Much of the theory developed for integral 
equations makes use of the £2 properties of the data and solution, but this is at 
odds with the system (2.8). As a higly non-uniform mesh is refined, the -^ 2 -iiorm of 
the vector a = (cri • • • cr^)^ becomes increasingly incomparable to the £2 norm of the 
solution to the continuous integral equation (2.7). For a set of quadrature weights {hj} 
which accurately integrate a and cr^, the proper discrete unknown should therefore 
be d-j = ^yhj (jj so that 


11 ^ 11^2 = 

3 

/ \ 2 

= Z] 

(2.13) 

3 

- j^<J^{x)ds{x) = ||cr||£ 2 - 


Intuitively, this embedding properly scales the unknown aj according to the clustering 
of the discretization along F. This re-weighting enables us to replace the discrete 
system in (2.8) with: 


(2.14) 


-\/hjCTj + Z 


\/hj Wje d 


Vhe dn^ 


■9k{xj,xe) ^/heae = -^/h~ju'’'^^{xj), 
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and declare & to be the new unknown. There is no reason to assume that the hj’s 
and the are the same, however in practice they are very similar except near the 

singularity of gk^ 

Under this re-weighting, the condition number of the discrete system converges to the 
condition number of the continuous problem as the mesh size tends to zero. If the 
curve r has corners, then, under refinement, the condition number of the original sys¬ 
tem (2.8) will usually diverge. For a thorough discussion and many results concerning 
this idea, see [7]. This norm-preserving embedding is one of the main tools used to 
construct high-order accurate boundary integral equation codes in complicated and 
singular geometries for both the Dirichlet and Neumann problems. Similar ideas with 
regard to £i-embedding have recently been used for divergence-form differential equa¬ 
tions with high-contrast background media [4]. Often these re-weighting techniques 
alleviate the need for designing specialized corner quadratures that are able to inte¬ 
grate singular Green’s functions multiplied by singular densities which diverge in the 
corner [12, 13, 11, 44]. 

It is with the previous section in mind that we begin to investigate the relationship 
between the solution of a scattering problem on a polygonal domain with that of 
a nearby smooth domain. In the next section we describe a simple convolution- 
based method to smooth polygons, and then report on the relationship between the 
numerical solutions to scattering problems in the smoothed and singular geometries. 

3. Smoothing polygons in 2-dimensions. An obvious approach to smooth¬ 
ing polygons is to locally represent the polygon as a graph and convolve with a smooth, 
compactly-supported even function with some specified order of differentiability. How¬ 
ever obvious, this technique seems not to have been analyzed or reported in the lit¬ 
erature. We use smooth to mean that the function is band-limited to some specified 
order. We restrict our attention to closed domains in two dimensions because of the 
emphasis on applications to scattering problems. Scattering from open surfaces re¬ 
quires several other numerical and analytical tools [15, 40, 41, 47]. Convolutional 
smoothing is an effective method in two dimensions due to the following elementary 
lemma: 

Lemma 3.1. Let (p(x) he an even, integrahle function, with compact support and total 
integral 1. For any a, 6 G M we have 


oo 


(3.1) 



Proof. This follows immediately from the observation that 


oo 


(3.2) 



— OO 


□ 

More importantly, this theorem remains true in n dimensions. If (/:? is a now a even 
function of n variables, with total integral 1, then a simple application of Fubini’s 
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(a) The rotated and translated graph of 
the neighborhood of a vertex of a poly¬ 
gon. 



0 ^ 


(b) The smoothed vertex of the graph 

fv- 


Fig. 1: The basic configuration for smoothing around a vertex. 


theorem shows that convolving (p with a linear function simply reproduces that func¬ 
tion. 

In what follows, let a polygon 7^ G be described by an ordered set of n +1 vertices 
{vj} and n edges {ej} such that vi = Each edge ej is defined by the set 

{vj, In a sufficiently small neighborhood of a particular vertex v, the polygon 

can be represented as an even graph of some function fy over a support line at v. We 
can normalize coordinates so that x = 0 corresponds to the vertex, with fy{0) = 0 . 
Then, for some (5 > 0, the function fy is linear on intervals [—(^,0] and [0,(5]. See 
Figure la for a plot of this configuration. Suppose that our convolution kernel (p is 
supported on [—1,1], then for some 0 < h < S/2 let 

(3.3) *>*(») dHD- 

and set 

h 

(3-4) f!^(x) = J (Ph{y) fv{x - y) dy. 

-h 


From the lemma, it is clear that 

(3.5) /’^{x) = fy{x) if \x\ > h. 

Hence the graph of defines a smooth (with band-limit dependent on that of p) 
curve that agrees with the graph of fy outside an neighborhood of the vertex of size 
h. See Figure lb for a depiction. 

This gives an effective means to smooth the vertices of the polygon V; since only a 
neighborhood of each vertex is changed, they can be smoothed locally and then glued 
together along the remaining straight edges. If the interior angle at a vertex is less 
than TT, then the smoothed vertex lies inside of the original polygon, whereas if it is 
larger than tt, then the smoothed vertex lies in the exterior. 
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Fig. 2: A range of smoothings of a 7r/2 corner done by convolving a local representation with '0^, 
with k = 8 andh = 0.025, 0.05, 0.1, 0.2, 0.4. 


The following simple algorithm can be used to uniformly smooth the polygon V with 
a given smooth, even function (^, with support in [—1,1]. 

Algorithm for polygonal smoothing via convolution 


Step 0: Choose a smoothing parameter h > 0, smaller 
than ^ min{|t;j — Vj-^i | : j = 1,..., n}. 

Step 1: For each represent a neighborhood of the vertex Vj as the 
graph of an even piecewise linear function fj over a support line 
to V at Vj. 

Step 2 : Convolve the functions fj with cph, to obtain f^. 

Step 3: Replace a neighborhood of Vj with part of the graph of by 
gluing along the linear parts of the graph of /j^, which agree 
with the graph of fj. 


Remark. The reason to use an even linear function in Step 1 is to insure that the 
smoothed polygon has the same discrete symmetries as V. 

The convolution can be done efficiently via either closed-form analytic expressions 
(depending on the choice of kernel cp) or by high-order numerical integration using an 
adaptive discretization scheme of the polygon and kernel as discussed in more detail 
in Section 3.3. Furthermore, an adaptive smoothing algorithm can be constructed 
by which the width parameter h is allowed to depend on the pairwise vertex spacing 
\Vj-Vj+l\. 


3.1. Selection of smoothing kernels. To make this an effective method re¬ 
quires the choice of a good family of smoothing kernels. We briefly discuss the details 
concerning two such kernels, one compactly supported and the other numerieally 
compactly supported. Let us first examine the family of functions f^k{x) G 

'^k{x) = Cfe (1 - X[-I,i\{x), 


(3.6) 
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(a) Plots of for ^ = 2, 4, 8, 16. 


- ^16 

(b) logj^Q of the absolute value of the Fourier 
transform of the kernels 



Fig. 3: Examples of the convolution kernels and their log-power spectra. 


where X[a,b] is the indicator function on the interval [a, 6]. These functions should be 
familiar from undergraduate analysis, and are well-suited to convolutional smoothing. 
Here is chosen so that ipk has total integral 1. In fact, 

(3^7) = + 0 Xif+i) 


An example of smoothing a right-angled vertex using this kernel is shown in Figure 2. 

When choosing a kernel with which to perform this convolutional smoothing, it is 
important to choose one which is localized in both physical space and Fourier space. 
Post-convolution, the resulting smooth curve will then have a band-limit proportional 
to the product of the band-limits of the straight edges and the kernel. The lower the 
resulting band-limit, the more accurately the curve can be discretized with a fixed 
number of degrees of freedom (discretization points). The Fourier transform of the 
function ipk is given analytically as 


(3.8) 


io = MO 


= r 




4+1 ( 20 , 


where Jn is the Bessel function of the first kind of order n and we have chosen the 
convention 


(3.9) 


/ OO 

fix) 

-C» 


It is clear that |'0/c(OI — '0/c(O) 


1, and asymptotically for large ^ these behave like: 


(3.10) 


|A(C)I 




This shows that once |^| > 2k/e, the Fourier transform of decays quite rapidly. 
The Fourier transform of the scaled function satisfies 





(0 =J^lA,h(x)j (C) = A(0> 


(3.11) 
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(a) Several smoothings of a triangle. 





(c) Several smoothings of a hexagon. 


(d) Several smoothings of a nonagon. 


Fig. 4: Convolutional smoothings of regular polygons. 


from which it follows that using frequencies a bit larger than 0{2k/eh) should suf¬ 
fice. Graphs of the Fourier transforms of {'^ 4 , '^ 12 , 'ipie} are shown in Figure 3. 

Figure 4 shows multiple smoothings of regular polygons convolved with the kernel 7 /;^ 
for various values of h. Note that the smoothings are nested inside one another for 
various values of h, with the more interior smoothings corresponding to larger values 
of h. 


The kernel ia equation (3.6) is convenient to use for our purposes because of its 
explicit compactness. However, if we are concerned with the support in the Fourier 
domain of (be. the band-limit of and therefore the band-limit of the smoothed 
geometry), we may wish to choose a kernel with somewhat more optimal uncertainty 
properties, the Gaussian: 


(3.12) 


0(6 =e-". 


The kernel (j) is not analytically compactly supported, however, it is numerically 
compactly supported. By this we mean that for any e > 0 we can find a threshold 
Xg > 0 such that for any \x\ > 0(x) < e. This, coupled with the integrability 
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Gaussian panels 



(a) The total geometry. 



Fig. 5: Smoothed polygon as sampled using Gauss-Legendre nodes. 


of 0 , allows us to choose a width parameter h such that outside of a neighborhood 
of a vertex, the resulting smoothed geometry differs pointwise from a straight line 
segment by at most e. Furthermore, if the neighborhood of a vertex is represented 
as the graph of a function /, the convolution of / with the Gaussian can be done 
analytically. Indeed, a symmetric / will be of the form f{x) = a\x\ + 6 , for some 
parameters a, 6 , and if we denote a scaled version of the Gaussian by (ph^ then 

= ax erf ( ) + 6 + 

\V2hJ 

where erf is the error function. Glearly, for any e > 0, there is a sufficiently large Xe 
such that |0/i * / — / I < e for all \x\ > x^. In the following numerical experiments, we 
set e ~ 10 “^^ such that smoothing calculations are done to nearly machine precision. 
It should be noted that the choice of e is independent of the choice of h. The value of 
e determines the size of | 0 /,,(h)|. 

3.2. Discretization of the smoothing. We first discretize a smoothed geom¬ 
etry with a specified value of h (depending on the particular polygon) using polyno¬ 
mial panels described by 16 Gauss-Legendre interpolation nodes (samples of values 
and derivatives are obtained numerically via adaptive discretization). Each panel is 
resolved when the corresponding Legendre polynomial coefficients (and those of the 
arclength function) of an oversampled discretization are below some threshold, set 
to 10“^^ in all cases. Obtaining higher precision is straightforward, and merely a 
matter of further refinement. We are mainly concerned with rough convergence on 
sub-wavelength rounded geometries. See Figure 5 for a picture of the discretization 
using Gauss-Legendre nodes on each panel, as well as a diagram of the smoothing 
kernel and corner. 

Outside of a distance h from the corner along an edge, the boundary contains straight 
edges which can be directly described using linear polynomial parameterizations. In¬ 
side a distance h from the corner, we insert (via translation and rotation) an adaptive 
panel-based discretization of the rounded function: 

(3.14) fg{x) = J (^a - dt, for X e {-w/2,w/2), 

where for e > 0, J = S{w) is chosen such that fs matches the original polygon to 
precision e. Figure 5 depicts the lengths a, re, and h. It is the curve fs that is 


(3.13) [0/.*/](x) 
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(a) The incoming field. 


(b) The scattered field. 


(c) The total field. 


Fig. 6: Example exterior sound-soft (Dirichlet) scattering problem. The real part of all fields is 
shown. The angle of incident plane wave is (p = —tt 14,. 


adaptively discretized so that its value, first derivative, and arclength functions are 
accurate to an absolute precision e [52]. In all examples, (j)s is the Gaussian kernel, 
and the explicit convolution is given in equation (3.13). In one final pre-processing 
step of the geometry, further refinement takes place until all neighboring panels differ 
in arclength by at most a factor of two and no panel is larger than 2A, where A 
is the wavelength inherent to the problem. Using the resulting discretization nodes 
we discretize the relevant integral equation (as in the next section) using the C 2 - 
weighted Nystrom method. This discretization scheme, used in conjunction with high- 
order quadratures for weakly-singular kernels, ensures the convergence of potentials 
for both the Dirichlet and Neumann scattering problems in corner geometries. 

We solve the linear system resulting from the Nystrom discretization of the continu¬ 
ous integral equation directly using the LAPACK implementation of LU -factorization. 
All numerical experiments are implemented in Fortran 90 and run using the Intel For¬ 
tran Compiler with MKL libraries. Entries in the discretized matrix corresponding 
to source-target pairs that reside on the same panel or on neighboring panels are 
determined using generalized Gaussian quadratures for logarithmically singular ker¬ 
nels [II]. Entries corresponding to source-target pairs that reside on non-neighboring 
panels are obtained from the 16-point Gaussian quadrature rule corresponding to unit 
weight (the Legendre polynomial case). 

3.3. Scattering from smoothed polygons: Sound-soft. We now turn our 
attention to numerical experiments pertaining to the scattering of acoustic waves from 
smoothed polygons. In this section, we study exterior Helmholtz scattering problems 
for Dirichlet boundary conditions; in the following section, we address the analogous 
Neumann problem. In the case of Dirichlet boundary conditions (corresponding to 
the case of a sound-soft scatterer), we have the following boundary value problem: 

=0 inM^V U, 

(3.15) ^ \ ^ 

^tot _ Q on F = 

along with suitable radiation conditions at infinity. Representing the scattered solu¬ 
tion u using a combined-field potential [26], 


(3.16) 


u = {Sk + i {ka + /3) Vk) cr, 
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we have the following second-kind integral equation along F for the density a: 

(3.17) I + (5* + z {ka + /3) %) a = -u*"" on T, 

where 5/c and are interpreted in their on-surface limiting sense. We have set a = 1.2 
and (3 = 0.8 in all examples. The scattered field is then calculated at all exterior 
volume locations using standard Gaussian quadrature for polynomials and the fast 
multipole method for the two-dimensional Helmholtz equation [30]. More accurate 
near-surface evaluation could be obtained using the methods of [35] or [43, 49]. 

The following simulations are obtained from driving the scattering problem by setting 
to be a two-dimensional plane-wave, traveling in the direction of the angle (j)\ 

(3.18) cos <l>+y sin 4>)_ 


It is easy to see that satisfies the free-space Helmholtz equation, but not the 
Sommerfeld radiation condition. See Figure 6 for depiction of an incoming plane wave 
scattered field and total field with Dirichlet boundary conditions. In this 
example, k = 12.43 + corresponding to a wavelength of A = 27r/ Rek ^ 0.505. 

The accuracy of the integral equation solver is tested by calculating the error in the 
potential when compared to a known solution obtained from placing a fundamental 
source in the interior of the object. Le., we solve a test problem: 


(3.19) 


(A + k‘^)u = 0 

u = gki’^xo) 


in \ 

on r. 


where Xq is placed near the center of the object. The potential u is then compared 
with the exact solution gk{'^Xo) at test points placed on a circle some distance away 
from the scatterer. 


We study the effect of the corner rounding by examining what is referred to as the 
sonar cross section (SCS) of the object O. Usually, this function is given in terms of 
the far-held behavior of the scattered held based on large-cc asymptotics of 


(3.20) «/“’-(a;) = ^ ^-ikr-y 

where r = x/\x\. The far-held signature is often used in inverse obstacle scatter¬ 
ing problems where measurement noise is frequently the dominant component any¬ 
way [21]. 

However, in our case, we have direct access to the scattered held at any observation 
point. We can thereby evaluate near-held functions at varying radii from the scatterer: 


= [ 9 k{d,y)a{y)ds{y), 

[6.21) Jr 

d = c-\- d cos Oi + d sin 0 j , 

where we denote the scattered held at a distance d from the centroid c of O. The 
vectors i, j are the unit vectors in the x, y directions, respectively. There are two types 
of cross sections that are usually computed: mono-static and bi-static. Mono-static 
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(a) The mono-static cross section. (b) The bi-static cross section. 


Fig. 7: Example mono-static and bi-static cross sections for the Dirichlet problem corresponding to 
the geometry in Figure 6, captured at approximately 20A from the origin. The absolute value of the 
scattered field is plotted on a log^^Q scale. 


cross sections characterize the scatterer in terms of the intensity of the backscatter 
in the same direction as the incoming wave. In particular, we calculate at a 

single value of 0 corresponding to the opposite angle of propagation of the incoming 
plane wave If the mono-static cross section is sampled at m angles, this requires 
solving m separate scattering problems. 

On the other hand, the bi-static cross section contains intensities of the scattered 
field for a fixed angle of incident plane wave. Figure 7 shows sample mono-static and 
bi-static cross sections for the scattering problem depicted in Figure 6, each captured 
at a distance of d = 10 ~ 20A from the origin. The angle of incidence for the bi-static 
case was 0 = —7r/4. In each case, the cross section is plotted on a polar grid in 
decibels: 

(3.22) C(0) = 10 logio (1^/(0)I). 

As the size of the region that is rounded near the corners is decreased, to below 
sub-wavelength, we see a convergence of the cross sections. Figure 8 shows a plot of 
several bi-static and mono-static cross sections for the same object (that in Figure 6). 
Here, we have increased the wave number to k = 54.32 + il0“^ to allow for a larger 
dynamic range of rounding widths. This value of k corresponds to a wavelength of 
A 0.12. The cross section is evaluated on a disc of radius 15 « 125A centered at the 
origin. 

The obvious question to ask is how close these solutions are to the solution in the 
case of scattering from an exact polygon with corners. Results of this experiment are 
shown in Figures 9, 10, 11, and 12. Convergence results of the far-held and moderately 
near-held bi-static cross sections are reported in Tables 10b, 11b, 12e and 12f. Near- 
held convergence is given in Table lid. In each case, the order of convergence of the 
scattered held is commensurate with the scale of the rounding. 

The errors in the value of the potential converge at a rate of roughly hrst-order with 
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(a) Several mono-static cross sections. 


(b) Details of several mono-static cross 
sections. 




(c) Several bi-static cross sections. (d) Details of several bi-static cross sec¬ 

tions. 


Fig. 8: Example mono-static and bi-static cross sections for the Dirichlet problem corresponding to 
several roundings of the geometry in Figure 6, captured on a disc of radius 15 ~ 125A from the origin. 
The absolute value of the scattered field is plotted on a decibel = lOlog^^Q scale. 


respect to the rounding parameter. Slightly faster convergence is actually observed, 
which may be due to the high accuracy of the rounding and the smoothing effects 
of the layer potential representation. We are currently investigating this phenomena. 
It is worth pointing out that in Figure 10 there are no correct digits in the solution 
(in a relative sense) until the rounding is performed on a scale roughly equal to the 
wavelength of the solution. The exact solution {h = 0.0) was calculated by dyadic 
refinement of the edges of the polygon near the corners to a scale of 10“^^. The 
resulting integral equation was solved using an £2 weighting scheme, as described 
in [7], 


3.4. Scattering from smoothed polygons: Sound-hard. We now present 
results corresponding to the sound-hard scattering problem, i.e. the exterior Neumann 
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(a) The empirical convergence for the comb 
shape with k = 54.32 + z 10“^ is and 

that for k = 6.79 + 1.2bi x IQ-® is 



(b) The empirical convergence for the triangle 
shape in Figure 12 with k = 52.13 + i is 
^(^1.19) ^ _ y. 77 + z 10“® is 


Fig. 9: Plot of the relative £2 error of the scattered potential for the Dirichlet problem versus rounding 
size in terms of wavelength. Both regimes exhibit commensurate convergence. 


10'^ 



(a) Errors in the real part of bi-static cross sec¬ 
tions. 
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(b) Errors in the bi-static cross sections. 


Fig. 10: Errors in the complex valued bi-static cross section for the Dirichlet problem at a distance 
of 1000 8333A from the origin (as compared to the true corner scattering problem). The error 

converges approximately to first order in the rounding parameter h. In each case, the PDF was 
solved to roughly a precision of 10“® in the Coo norm (as determined by testing against a known 
solution). 
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(a) Total field for k = 6.79 + 1.25z x 10“ 
wavelength is approximately A .93. 


\ The 
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(b) Errors in the bi-static cross sections for 
k = 6.79 + 1.25i x at a radius of 1000 
1081A. 



(c) Total field for k = 54.32 -\-i x 10“^ along with 
testing curve for convergence of scattered field. 
The corners in this plot were rounded with a pa¬ 
rameter h = 0.025. 
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(d) Errors in the bi-static cross section for 
k = 54.32 + z X 10“^ in the near-held at a 
radius of 3.5 30A. 


Fig. 11: Errors in the complex bi-static cross section for the Dirichlet problem (as compared to the 
true corner scattering problem). The error converges approximately to hrst order in the rounding 
parameter h. In each case, the PDE was solved to roughly a precision of 10“^ in the Cog norm (as 
determined by testing against a known solution). 





































































h 

n 

RMSE error 

Rel. ^2 error 

0.4 

1056 

5.0 X 10 “^ 

2.2 X 10“i 

0.2 

1248 

2.4 X 10-2 

1.1 X 10-1 

0.1 

1344 

1.1 X 10-2 

4.6 X 10-2 

0.05 

1440 

4.5 X 10 -^ 

2.0 X 10-2 

0.025 

1632 

1.9 X 10 -^ 

8.5 X 10-3 

0.0125 

0.0 

1728 

3456 

8.4 X lO -'^ 

3.7 X 10-3 


(e) Errors in the bi-static cross section at r 
10 for k = 52.13+ a0-'^. 
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(f) Errors in the bi-static cross section at r 
10 for k = 7.77+ 


Fig. 12: A depiction of sound-soft (Dirichlet) scattering for various roundings, along with convergence 
of the bi-static cross section in the moderate near-field. In each case, the PDE was solved to roughly 
a precision of 10“^ in the Coo norm (as determined by testing against a known solution). 
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(a) The incoming field. 


(b) The scattered field. 


(c) The total field. 


Fig. 13: Example exterior sound-hard (Neumann) scattering problem for k = 12.43 + zlO ^. The 
real part of all fields are shown. The angle of the incident plane wave is (/) = —7r/4. 


problem for the Helmholtz equation: 


(3.23) 


(A + = 0 

in \ fl. 


dn 

on r = on 


along with suitable radiation conditions at infinity. Representing the scattered solu¬ 
tion u using a single-layer potential: 

(3.24) u = Sk<T, 

we have the following second-kind integral equation along F for the density a: 

(7 

(3.25) __+5;a = -— onr, 

where = dSk/dn and is interpreted suitably as an on-surface limit. As before, our 
reference solver for the true corner problem follows the method detailed in [7]. 

We also recall that using representation (3.24) may yield spurious resonance in the 
resulting integral equation for values of k which correspond to eigenvalues of the 
interior Laplace Dirichlet problem. For simplicity we have chosen k to avoid these 
values. Well-conditioned combined-field representations exist which are invertible for 
all values of k with Im/c > 0, but they involve the composition of layer potentials, 
as in (2.9), not merely the summation [21]. After solving (3.25), we evaluate the 
scattered field as in the previous section, using the fast multipole method for the 
two-dimensional Helmholtz equation and standard Gaussian quadrature. 

The following simulations are obtained from driving the scattering problem by setting 
to be a two-dimensional plane-wave, as before, traveling in the direction of the 
angle 0: 

(3.26) uf‘={x) = 

See Figure 13 for a depiction of an incoming plane wave scattered field and 

total field u^^^ with Neumann boundary conditions. In this example, k = 12.43+H0“^, 
corresponding to a wavelength of A = 27r/ Re /c ~ 0.505. 
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The accuracy of the integral equation solver is tested, as before in (3.19), by com¬ 
parison with a known test solution. In order to study the effect of corner rounding 
for the Neumann problem, we reproduce several of the experiments performed in the 
Dirichlet case. In particular, we compare the bi-static SCS of the true corner problem 
with that from successive roundings. See Figures 14, 15, and 16 for plots of Neumann 
solutions and convergence results. 

As in the Dirichlet case, as the size of the region that is rounded near the corners is 
decreased, to below sub-wavelength, we see a convergence of the bi-static cross section 
of roughly first-order. We simulated the Neumann problem at the same frequencies 
as in the Dirichlet case for comparison. 

3.5. Extension to piecewise smooth boundaries. This technique can also 
be extended to piecewise smooth curvilinear polygons. Since we need a variant of this 
idea to smooth polyhedra in we pause to briefly describe it here. In short, in the 
neighborhood near a geometric singularity it is possible to construct a diffeomorphism 
to a truncated cone. The corner rounding can then be performed on the polygonal 
cone, and finally composed with the inverse of the diffeomorphism to smooth the 
original curvilinear polygon. 

To this end, let 7^ be a region in M whose boundary is composed of a finite collection 
of smoothly embedded arcs, {yj : j = 1 ,... ,n} meeting at points 

(3.27) Vj = n 7j+i 

and angles {0 < Oj < 27r}. We let 7 n+i denote a second copy of 71 . We are excluding 
the case of a cusp, i.e. Oj = 27r. 

Once again the idea is to change only a small neighborhood each vertex. We define 
(in complex notation) the planar regions 

(3 28) Wj = {z:Q<avgz<ej anA\z\<l}, if < tt, 

Wj = {^ : 0 < arg z <2tt — Oj and |^| < 1 }, if Oj > tt. 

Suppose that for each j for which Oj < tt we can find a diffeomorphism ipj from Wj 
to a neighborhood of Vj in 7 ^, which carries: 

0 —y Vj , 

(3.29) {z : arg^ = 0} D dWj a ray in 7 ^, 

{z : arg z = 0j}C^ dWj a ray in 7 ^+ 1 . 

> TT, then 'i^j is defined from a neighborhood of the vertex in Wj to a neighborhood 
of Vj in 7 ^^, with the boundary correspondence as before. Conformal mapping provides 
one effective method to define such maps. Other, more elementary techniques are 
also available. One such method, which works for regions with convex boundaries, is 
described in Section 6 . 

For each h > 0 we define W^ as the regions obtain by smoothing the vertex of Wj at 0 
as described above. For small h, we have W^ C Wj^ and the boundaries of W^ and 
Wj coincide outside of a small neighborhood of 0. Thus, for small enough r, the image 
^ljj{dW^ \ 5r(0)) lies along the boundary of V outside a small neighborhood of Vj. 
Hence the image ipj{dW^) defines a smoothing of the vertex at Vj. This procedure 
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(b) Errors in the bi-static cross sections for 
k = 6.79 + 1.25i x at a radius of 1000 
1081A. 


(a) Total field for k = 6.79 + 1.25z x lO”®. The 
wavelength is approximately A .93, and the 
rounding parameter was h = 0.2. 
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(d) Errors in the bi-static cross section for 
k = 54.32 + z X 10“^ in the near-field at a 
radius of 3.5 ~ 30A. 


(c) Total field for k = 54.32 + z x 10“^ along with 
testing curve for convergence of scattered field. 
The corners in this plot were rounded with a pa¬ 
rameter h = 0.025. 


Fig. 14: Errors in the complex bi-static cross section (as compared to the true corner scattering 
problem) for the Neumann problem. The error converges approximately to first order in the rounding 
parameter h. In each case, the PDE was solved to roughly a precision of 10“® in the Coo norm (as 
determined by testing against a known solution). 
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(c) Errors in the bi-static cross section at r 
10 for k = 52.13+ a0-'^. 
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(d) Errors in the bi-static cross section at r 
10 for k = 7.77+ 


Fig. 15: A depiction of sound-hard (Neumann) scattering for various roundings, along with conver¬ 
gence of the bi-static cross section in the moderate near-field. In each case, the PDE was solved to 
roughly a precision of 10“® in the £00 norm (as determined by testing against a known solution). 



(a) The empirical convergence for the comb 
shape with k = 54.32 + z 10“^ is and 

that for k = 6.79 + 1.25z x 10“® is 



(b) The empirical convergence for the triangle 
shape in Figure 15 with k = 52.13 + i 10~^ is 
and that for k = 7.77 + z 10“® is 

(!)(hi-20). 


Fig. 16: Plot of the relative £2 error of the scattered potential versus rounding size in terms of 
wavelength for the sound-hard (Neumann) problem. Both examples exhibit similar convergence. 
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is done locally in a small neighborhood of each vertex, allowing one to smooth the 
vertices while leaving as much of the remainder of the boundary of V fixed as desired. 


4. Polyhedra in three dimensions. In this section we describe several meth¬ 
ods for smoothing piecewise smooth boundaries of regions in In Section 4.1 we 
describe a special class of polyhedra, 3-regular Hamiltonian polyhedra^ whose bound¬ 
aries can be smoothed using the method described above with a parameter. In fact, 
all convex polyhedra, and many non-convex polyhedra can be smoothed this way, 
but the results are often not-optimal. In Section 4.2 we show that by modifying a 
polyhedron in a small neighborhood of its vertices one can obtain a 3-regular, Hamil¬ 
tonian polyhedron. Hence it can be smoothed using the method given in Section 4.1. 
This leads to a smoothed boundary that agrees with the original polyhedron outside 
a small neighborhood of the original edges and vertices. Unfortunately, the smoothed 
polyhedron will also contain open subsets of translated support planes of the vertices. 
This is both unsightly and can produced a dramatically enhanced scattered wave in 
the direction normal to the plane. A more robust approach is described in Section 4.3. 

4.1. 3-Regular Hamiltonian Polyhedra. There is a special collection of poly¬ 
hedra in whose edges and vertices can be smoothed using only what might be called 
the two-dimensional method with parameter. We first define this class: 

Definition 4.1. Let P he a polyhedron in and Gp the graph defined by its edges. 
P is 3-regular if every vertex is the interseetion of three faees, or, equivalently, ifGp 
is a 3-regular graph. It is Hamiltonian if there is a finite eolleetion of disjoint eyeles 
{Cl ,... ,Gi} (Z Gp so that every vertex belongs to exaetly one of these eyeles. 

It turns out that not every 3-regular polyhedron is Hamiltonian, and when one is, 
the problem of finding these cycles is not generally solvable in polynomial time. On 
the other hand, all 3-regular Platonic solids (tetrahedron, cube, and dodecahedron) 
are Hamiltonian, as well as many examples that arise in practice. As this class of 
polyhedra can be smoothed by smoothing only edges, we take a moment to describe 
the procedure. 

Let P be a 3-regular Hamiltonian polyhedron, with C = {Ci,..., C/} a collection of 
disjoint cycles exhausting the vertices. Let £ = {ei,..., e^} be the edges of Gp that 
are not contained in any cycle. Because the graph is 3-regular, we know that every 
vertex in Gp lies on exactly one of these edges. Moreover the edges in £ are disjoint. 

To smooth P we first smooth the corners that lie along the edges in £ (using the 
two-dimensional method described earlier). It is easy to see that every edge ej lies 
in the intersection of two planes {7Tj^,7rj^}. Let be a plane orthogonal to the line 
£j = TTj^ n TTj^, and ujj = {cp, 0) G S‘^ be the direction of £j. Finally let Zj denote the 
component of (tt^^ U ) D Tr^g so that a neighborhood of ej in P lies inside the corner 


(4.1) 


Kj = {qP • I ^ Ij t G M}. 


If 7 ' is a smoothing of this curve, as defined in the previous section, then 


(4.2) 


N = {9 + i e R} 
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is a smoothing of the corner. If the rounding of 7 ^ is done close enough to the vertex, 
then we can smoothly replace a neighborhood of Cj in P with its smoothed version 
in Kj by simply intersecting the interior of the region bounded by Kj with P. Away 
from the smoothed edge, Kj is still a union of planar regions which can be glued onto 
P, thereby replacing ej with a smooth transition between these planar regions. 

Since the edges in £ do not intersect, each of these smoothing operations can be done 
independently of the others. Let P' denote the body in obtained by smoothing all 
of these edges. Since every vertex lies on one of the edges in f, the cycles on P are 
replaced by cycles C' = {C[^ ..., C[} on P' that are smooth non-intersecting curves. 
That is to say, every vertex has been smoothed. The boundary of the body P' is 
a comprised of bounded smooth surfaces, which are mostly planar regions. These 
surfaces are bounded by smooth, disjoint, closed curves, along which these surfaces 
meet. All that remains is to smooth these curves of intersection. 

To that end we now define a diffeomorphism from a standard model onto a neighbor¬ 
hood of Cj. We smooth the standard model and use this map to glue the result into 
P'. Let Cj : [0, Lj] P' be an arclength parameterization of Cj. The unit vector field 
Tj(t) = dtCj is tangent to Cj. Two smooth surfaces Sij and S 2 j meet, transversely, 
along this curve. Let Nij(t) be the unit vector normal to Tj(t) lying along Sij^ i = 1 , 2 . 
Let 7rj{t) denote the plane through Cj{t) spanned by {Nij{t), N 2 j{t)}. 

For e > 0, let Uje denote the e-neighborhood of Cj. There is a radius e > 0 so 
that these planes {7rj{t) : t G [0,Pj]} define a foliation of Uje. This follows from the 
inverse function theorem and the compactness of the curve. We define a map from 
Vjs = [0,I/j) X {—S,6) X {—S,6) into a neighborhood of Cj by letting 

(4.3) ^j{t, 5i, 52 ) = Cj{t) + siNij{t) + S2N2j{t). 

The differential of at (t, 0, 0) is given by 

(4.4) (i<Fj(t,0,0) = Tj{t)dt P Nij{t)dsi P N2j{t)ds2, 

which is clearly of rank three. From the inverse function theorem it now follows easily 
that there is an ^ > 0 so that \y.^ is a diffeomorphism onto its image, which is a 
neighborhood of Cj foliated by the planes {'T^j{t)}. We can continue as a smooth 
Pj-periodic function. 

We first make the assumption that P' fl Uje lies in the image of the positive orthant 
in the (si, 52 )-variables under this map. This is certainly the case if the interior angle 
along Cj is everywhere less than tt. Under this assumption it is easy to see that 
for an 7 > 0 there is a set of the form Vjr^ = M x [ 0 ,??] X [0, 77 ] on which is a 
penodic-diffeomorphism. Moreover 4>j(I^-^) C P' exhausts a neighborhood of Cj in 
P'. 

We now smooth the corner Vjj^ to obtain where the smoothed edge lies in . The 
image of the smoothed corner under ^j defines a smoothing of Cj. As these curves 
are disjoint, each one can be smoothed independently. We let P" denote the resulting 
body in It is a smoothed approximation to P. If the interior angle is greater than 
TT at every point, then we can smooth the corner by smoothing the exterior, which 
satisfies the hypotheses above. 

We demonstrate this approach to smoothing polyhedra by smoothing a cubic torus P. 
We suppose that P is oriented parallel to the standard coordinate axes. There are four 
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(a) Rounded edges and cycles of the cubic (b) Rounded edges/cycles of a face, 

torus. 


Fig. 17: Rounding Hamiltonian cycles of a cubic torus. Images were constructed merely for illustra¬ 
tive purposes only using FreeCAD. 


cycles {Cl,..., C 4 }, each containing four edges and parallel to the xz-plane. These 
cycles bound the faces that have non-trivial topology. For the edges not belonging to 
cycles, f, we use the eight edges parallel to the ^-axis. If the edges in £ are smoothed, 
then the cross sections of P' perpendicular to the ^-axis are smoothed squares, as 
shown in Figure 17. 

We now smooth the remaining edges using the representation in equation (4.3) along 
with the smoothing of the right angle used to smooth the edges in £. Figure 18 
shows two views of the upper part of the final smoothed cubic torus. We should point 
out that these images were constructed using the software FreeCAD using low-order 
fillet procedures, and are illustrative only. Constructing the high-order computational 
geometry software to carry out convolutional smoothing and subsequent high-order 
piecewise triangulation for polyhedra in three dimensions is an ongoing project. 

4.2. Smoothing the Vertices: I. The method for smoothing edges described 
in the previous section can be used to smooth an arbitrary convex polyhedron in a 
two-step procedure. Let 7^ be a convex polyhedron with faces T — {/i,..., //}, edges 
£ = {ei,..., e^}, and vertices V = {i;i,..., At each vertex we choose an outward 
pointing support vector, {i/i,..., i/^}. Suppose that the edges at the vertex join 
to the vertices (t’/ci, • • •, A good choice for Uj is to take 


1 P 



(4.5) 


as it will preserve whatever symmetries the original polyhedron possesses in the 
smoothed domain. 

Given e > 0 , we define a neighborhood of the vertices by the condition 


(4.6) 


X G K if for some j we have (X — vj, tyj) > —e. 
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(a) Final smooth surface. (b) Smoothed exterior corner. 


Fig. 18: Final rounding of Hamiltonian cycles of a cubic torus. Images were constructed merely for 
illustrative purposes only using FreeCAD. 



(c) Smoothed interior corner. 


Note that V C Vq, and, for small e > 0, the intersection dV^ fl 7^ is a disjoint union 
of small polygons lying near the vertices. See Figure 19a. It is easy to see that 
the resultant polyhedron is 3-regular and Hamiltonian, with the disjoint cycles being 
those introduced when cutting off the vertices. 

To smooth the polyhedron we first smooth the edges. An edge lies in the inter¬ 
section of two faces D . Let iTk be the plane, through the midpoint of the edge, 
which is perpendicular to e/^. Using the method described in Section 3 we can smooth 
the vertex Htt/c of the polygon defined by the intersection of tt/^ with V. By parallel 
translating this smoothed vertex along the edge, we can replace a neighborhood of 
the edge by a smooth surface joining the plane containing fi^ to the plane contain¬ 
ing With h > 0 the smoothing parameter from Section 3, we let Vh denote the 
polygon with all its edges smoothed in this manner. 

Of course, near enough to a vertex, the smoothings of different edges intersect, but 
given e > 0 we can choose a sufficiently small h > 0 so that, in the set U/, the 
modifications corresponding to the different edges are disjoint. With such choices, 
the intersection dV^ fl 7^/^ is a disjoint union of polygons with smoothed vertices lying 
in the planes 

(4.7) {X=-€. 


Using the technique described in Section 4.1 the edges along which these smoothed 
polygons meet dVh can be smoothed, leading to an overall smoothing of the original 
polyhedron. While it is clear that this can be done in an arbitrarily small neigh¬ 
borhood of the singular locus of dV, the vertex is replaced by a smooth surface 
containing a open subset of the plane defined in (4.7). For applications to scattering 
theory this might not be desirable, as it will produce a considerable amplification of 
the scattered signal in this direction. In the next section we describe a method for 
smoothing vertices that produces a better result. 

4.3. Smoothing the Vertices, II. The second method for smoothing vertices 
takes as its starting point the domain Vh constructed in the previous sections by 
smoothing the edges (as depicted in Figure 19b). We assume that there is a positive 
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(a) The intersections of V with planes 
{X-Vj,Uj) = -e. 



(b) View along the 2 ;-axis of the intersections 
of the smoothed edges with planes 
ix = -e. 


Fig. 19: Corner rounding of a tetrahedron. Images are for illustrative purposes only, constructed 
using FreeCAD. 


eo SO that the intersections 

(4.8) = p;, n {X : (X - = -eo} 

are disjoint smoothed polygons. With this assumption each vertex can be smoothed 
without reference to any other vertex. We can therefore fix a j and describe the 
method for smoothing 7^ in a neighborhood of v^. 

We let 0 < ej < eo denote the infimum of the numbers so that is a polygon with 
smoothed vertices. The domain is a smoothed polygon, where the smoothings 
of two (or more) of the edges meet without any flat segment between them. 

For each e > e^, we let denote a maximally smooth parameterization of dP^^^ on 
the unit circle. That is, is a map from to dP^^P Therefore, we can represent 
it in terms of a Fourier expansion: 

CX3 

(4.9) i>e(0)= ^ 

n= — oo 

The infinite sum defines a map from the unit circle to dP^^^ translated to the plane 
(X, Uj) = 0. We adjust eo so that eo > 4ej for all j. 

For each Cj < e < cq we can extend this map as a diffeomorphism from the unit disk 
to the smoothed polygon P^'P For example, since the boundary of P^'^ is convex, 
it follows from a theorem of Choquet that the harmonic extension has the desired 
properties: 

c» 

(4.10) $^(r,6>)= - ej/j. 

n=—oo 

For additional details, see the next section and [20]. The image of 4>e lies in the 
plane {X — Vj^ Vj) = —e. 
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To use these maps to define a smoothing we need to choose two auxiliary functions. 
First we choose a number rji so that 2ej < rji < eg. Next, choose a smooth, convex, 
increasing function x(s) defined in [ 0 , eg] so that for s > 771 , x(s) = s, x( 0 ) = e^, and 

(4.11) xf"^^(0)=0 for m = 1,...,/c. 

We also choose positive numbers rg < eg, 772 < '^o/2, and an even convex function 
'ip{r) defined in a neighborhood of 0. We require 

'0(r) = r for r > rg, 

(4.12) 'ip{0)=r]2, 

^[m] _ Q for 772 = 1, . . . , /c. 

The smoothing of the neighborhood is defined as the image of [0, eg] x under the 
map 

(4.13) 4'(r, 0) : (r, 0) ^ “ x{r)i'j- 

For r > max{77i,rg} this map simplifies to 

(4.14) (r, 0) i-A (1, 0) — ruj. 

That is to say, its image lies in the already smoothed part of dVh near to Vj. Our 
assumptions assure, that as function of x = rcosO and y = rsin^, the map {x^y) i-A 
T(x,^) is at least 772 — 1 times differentiable in a neighborhood of (0,0) and dT(0,0) 
has rank 2. Therefore, the image of ^^^(O) under T is a smooth sub-manifold of M^. 
The image lies in the set 

(4.15) {X - <-Cj, 

where the vector 12 j is the normal vector to the smoothed vertex at the point T(0, 0). 

In describing this method for smoothing vertices, we have assumed that the original 
polyhedron is convex, but this is not necessary for the method to be applicable. It is 
merely required that each vertex Vj has a local strict supporting plane. This means 
that there is a vector Uj so that if (uj, Uj) = Cj, then for some r > 0, 

(4.16) r n Br{vj) \ Vj C {X : {X, Uj) < Cj}. 

Here Br{v) = {X G : |X — v| < r}. The existence of a strict local support plane 
implies that for a range of e > 0 the sets 

(4.17) Pi =Vn{X:{X-Vj,i.j)=-e} 

are polygons. With this assumption we proceed as before, first smoothing the edges 
to produce Vh- For some 0 < eg < e < ei the sets 

(4.18) =rhn{X:{X- Vj,Vj) = -e} 

are smoothed polygons. The method described above can easily be adapted to smooth 
the vertex in this case as well. The results of using this technique to smooth polyhedra 
are shown in Figure 20. 




Fig. 20: Several smoothed polyhedra. 


5. General Polyhedra. Using this general scheme of first smoothing the edges, 
and then using diffeomorphisms to smooth the vertices we now describe a method 
that suffices to smooth arbitrary globally embedded polyhedra in Let P be a 
polyhedron, by which we mean a bounded region in whose boundary is a union of 
polygons lying in planes. We let {vj} denote the vertices of V. If the polyhedron has a 
strict local support plane at every vertex, then the method describe in Section 4.3 can 
be applied to produce a locally smoothed polyhedron, by first smoothing the edges 
and then the vertices. There are polyhedra that do not have strict local support 
planes at every vertex, e.g. the cubic torus does not have support planes at the inner 
vertices. 

The method of Section 4.3 requires that near to the vertex Vj^ the polyhedron is the 
cone over an intersection with a plane of the form 

(5.1) Vr\{X-.{X-v^,u^) = -e}. 

Let Sr{v) denote the sphere of radius r centered at v. A slightly more complicated 
method results if we instead assume that for each j there is an rj > 0 so that 

1. P n Srj{vj) is a connected region Rj on Srj{vj) bounded by a simple closed 
curve, 7 j, 

2 . V n Brj{vj) is the cone over Rj with vertex Vj. 

The curve 7 j is a piecewise geodesic polygon on the sphere. A polyhedron satisfying 
these conditions is globally embedded. In general the region Rj could have several 
connected components, a case that we do not consider further. 

For sufficiently small h > 0 , we let Vh denote the result of smoothing the edges of V 
as described in Section 4.2. If V is globally embedded, then, at each vertex there is a 
range of radii poj{h) < r < pij{h) so that the intersections 

(5.2) VhnSrivj) 

are regions Rj{r,h) bounded by simple closed curves, jj{r^h), that are smooth¬ 
ings of the curves V fl Sr{vj). As h ^ 0 , it is clear that poj{h) tends to 0 and 

\iminfh^o Pij{h) > rj. 






C. L. EPSTEIN AND M. O’NEIL 


31 


For each poj{h) < r < pij{h) we let 

(5.3) : Di{0) —^ Rj{r,h) 

be a diffeomorphism from the unit disk onto the region Rj{r, h). The maps {^r} can, 
for example, be defined as the conformal maps from T)i(O) to the spherical domain 
Rj{r, h), normalized so that 0 is mapped to points lying on a carefully selected curve. 
Using these maps we can define an analogue of the map T(r, 6>), defined in (4.13), so 
that the image of [0, eo] x under this map is a smoothed version of a neighborhood 
of the vertex vj, which is joined smoothly to Vh- We leave the detailed construction 
of these maps to the ambitious reader. 

There are also approaches to smoothing that first smooth the vertices, using the 
methods described above, and then interpolate these smoothings along the edges. It 
is very difficult to preserve convexity using this order of operations. That is why we 
have only described methods that first smooth the edges, and then the vertices, using 
a slicing approach along with families of diffeomorphisms. 

This completes the description of our algorithms for smoothing polyhedra in Note 
that one can restrict the modifications of the original polyhedron to lie in an arbitrarily 
specified neighborhood of the 1-skeleton of the boundary of P. One also retains consid¬ 
erable control on the relationship between the Gauss map of the smoothed polyhedron 
and that of the original, which is crucial for the behavior of scattered waves. In the 
final section we provide several practical methods for constructing diffeomorphisms. 

6. Methods to construct diffeomorphisms. We now describe several meth¬ 
ods to define extensions of a map from to T, a Jordan curve in the plane, which are 
diffeomorphisms from the unit disk Pi(0) to the region Dr which is bounded by T. 

6.1. Method 1. Conformal mapping provides a method that can be computa¬ 
tionally expensive and numerically ill-conditioned (depending on the geometry) [46], 
but guaranteed to work in considerable generality. In particular Dr can be a simply 
connected region in either a plane or a round sphere. Suppose that / : Pi(0) —> Dr 
is a conformal diffeomorphism. If T is convex, then 

(6.1) Tr = {fire^‘>): 9G[0,2n]} 

is a convex curve for every r G (0,1]. If T is star shaped with respect to 0 and / is 
normalized so that /(O) = 0 then the curves {T^ : r G (0,1]} are star shaped. These 
results can be found in [48]. 

6.2. Method 2. There is a simple method that is guaranteed to give a diffeo¬ 
morphism if Dr lies in a plane and the initial curve T is convex. A theorem of T. 
Rado, Kneser, and G. Choquet states that if (i4, v) defines a homeomorphism from 
the unit circle to T which bounds a convex region Dr, then the harmonic extension 
of the coordinate functions {U,V) defines a diffeomorphism from the interior of Di 
to Dr [20]. This theorem does not require T to be strictly convex or smooth. 

If the boundary map is given in terms of the Fourier series 

( CX3 OO \ 

j=-oo j=-oo J 


(6.2) 
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it follows from Choquet’s theorem that 


( oo oo \ 

j=-oo j=-oo J 

defines a diffeomorphism from Di{0) onto I^r- 

6.3. Method 3. If we specify a convex curve F in terms of its Gauss map, that 
is, as the image 

(6.4) G{0) = g{0) {cos 0, sin 0) g'{0) (—sin6>, cos6>), 

then we can proceed as above to get a diffeomorphism. If g has the Fourier represen¬ 
tation 

oo 

(6.5) g{e)= /3ne'"^ 

n= —OO 

then we can again apply Choquet’s theorem to construct a harmonic extension, which 
is guaranteed to give a diffeomorphism. The map defined in (6.4) can be represented 
as 

(6.6) e*® n- {g{0) + ig'{8))e'-^. 

Using the Fourier representation in (6.5), we see that 

oo 

(6.7) G{d,r)= ^n-i(2-n)rHe-^ 

n= —OO 

defines the harmonic extension of this map, and is therefore a diffeomorphism from 
1^1 (0) onto I^r- 

7. Conclusions. In this paper we have presented several algorithms for modi¬ 
fying polygons and polyhedra into fully regularized surfaces without geometric sin¬ 
gularities (vertices and edges). The original polygon (or polyhedron) is modified in a 
controllable and arbitrarily small neighborhood of its singular set. We have compared 
the solution to acoustic scattering problems from the original singular boundary to 
that obtained by smoothing the boundary at sub-wavelength scales in two dimensions. 
Both near- and far-held solutions converge at a rate slightly faster than hrst-order in 
the rounding parameter. Understanding this rate of convergence is an ongoing re¬ 
search topic in our group. 

Constructing numerical codes for performing rounding in two dimensions is relatively 
straightforward. We presented results for the polygonal case; software implementing 
the rounding of vertices joining piecewise smooth curves is currently under develop¬ 
ment, requiring merely the re-parameterization of the curve near the vertex as a graph 
above a support line tangent to the vertex. These computations are relatively fast, 
efficient, and accurate to near machine precision in two dimensions. 

We also introduced the analytical foundation for constructing high-order roundings 
of polyhedra in three dimensions. Composing various methods with diffeomorphisms 
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near vertices allows for similar regularizations to be computed as in the two-dimensional 
case. Building more efficient software to perform these computations is a work in 
progress. 

Preliminary Matlab code which performs the vertex and edge smoothing for convex 
polyhedra in three dimensions has been made available at: 

http://gitlab.com/oneilm/rounding 

If only approximate scattering solutions are required to the true problem involving 
geometries with corners and edges, the algorithms of this paper offers a method to 
obtain these results with reduced computational cost and controlled accuracy. Fur¬ 
thermore, the methods require nothing other than the usual quadratures for weakly- 
singular functions on smooth curves or surfaces. Full extensions of these smoothing 
algorithms to three dimensions may have a wide array of applications in high-order 
CAD and CAE packages, as many existing software solutions only allow for twice 
differentiable roundings (fillets). 

Lastly, we would like to note that the algorithms presented in this paper for geomet¬ 
ric regularization in three dimensions are only one piece of a larger effort to develop 
high-order scattering codes for arbitrary geometries. In three dimensions, all the 
numerical tools that are required to solve boundary integral equations are more ex¬ 
pensive and more sophisticated than those in two dimensions. Merely constructing 
high-order Nystrom-compatible quadratures for the function l/|cc — y\ along trian¬ 
gular patches is a relatively recent result [9, 10]. Coupling these schemes with fast 
algorithms and high-order triangulations is under active development. Performing the 
analogous convergence studies for rounding in three dimensions will be reported at 
a later date, after the requisite high-order accurate computational PDE algorithms 
have been developed. 
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